Quantum Computation of a Complex System: the Kicked Harper Model 
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The simulation of complex quantum systems on a quantum computer is studied, taking the kicked 
Harper model as an example. This well-studied system has a rich variety of dynamical behavior de- 
pending on parameters, displays interesting phenomena such as fractal spectra, mixed phase space, 
dynamical localization, anomalous diffusion, or partial delocalization, and can describe electrons 
in a magnetic field. Three different quantum algorithms are presented and analyzed, enabling to 
simulate efficiently the evolution operator of this system with different precision using difi'erent re- 
sources. Depending on the parameters chosen, the system is near-integrable, localized, or partially 
delocalized. In each case we identify transport or spectral quantities which can be obtained more 
efficiently on a quantum computer than on a classical one. In most polynomial gain com- 

pared to classical algorithms is obtained, which can be quadratic or less depending on the parameter 
regime. We also present the effects of static imperfections on the quantities selected, and show that 
depending on the regime of parameters, very difi'erent behaviors are observed. Some quantities can 
be obtained reliably with moderate levels of imperfection, whereas others are exponentially sensi- 
tive to imperfection strength. In particular, the imperfection threshold for delocalization becomes 
exponentially small in the partially delocalized regime. Our results show that interesting behavior 
can be observed with as little as 7-8 qubits, and can be reliably measured in presence of moderate 
levels of internal imperfections. 
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I. INTRODUCTION 



In the past few years, the field of quantum information 
has attracted more and more attention in the scien- 
tific community. Among the most fascinating promises 
of this domain is the possibility of building a quantum 
computer. Such a quantum processor can use the su- 
perposition principle and the interferences of quantum 
mechanics to perform new types of algorithms which can 
be much more efficient than classical algorithms. Cele- 
brated examples are Shor's algorithm which factors large 
integers exponentially faster than any known classical al- 
gorithm j3| and Grover's algorithm which searches un- 
structured lists quadratically faster than classical meth- 
ods Another type of quantum algorithms concerns 
the simulation of physical systems. Examples include 
qua,] 

spin systems |3 , 
rithms implementing quantum maps are especially inter- 
esting, since the systems simulated have simple equations 
of motion but can display very complex behaviors. Their 
simplicity enables to simulate them with a small number 
of qubits. For example, it is possible to simulate effi- 
ciently the baker map j8| (experimental implementation 
with the NMR technique has already been performed |^), 
the quantum kicked rotator , the sawtooth map 

p^, or the tent map T^. In such algorithms, it is im- 
portant to determine which physical quantities can be 
obtained accurately through measurement on the quan- 
tum computer, and what is the total algorithmic com- 
plexity of the whole process. It is equally important to 
determine the effects of errors in the computation to as- 
sess the efficiency of the algorithm on a realistic quantum 



many-body quantum systems ,41 , classical and quantum 
I , classical dynamical systems |^ Q . Algo- 



computer. 

In the present paper, we will study in detail an im- 
portant example of quantum map, namely the kicked 
Harper model. The Hamiltonian of this system has a 
simple form, yet displays many interesting physical fea- 
tures not present in quantum maps previously studied 
in this context, such as fractal spectra, stochastic web, 
anomalous diffusion, or coexistence of localized and de- 
localized states. It was introduced in the context of solid 
state physics (motion of electrons in presence of mag- 
netic field), and has been the subject of many studies. 
Using this model as a test ground, we will present three 
different ways of simulating the quantum map on a quan- 
tum computer, two of them inspired by previous works, 
and compare their efficiency. We will then present ex- 
amples of physical quantities which can be obtained on 
a quantum computer. It turns out that depending on 
the parameters of the system, at least polynomial speed- 
up compared to classical algorithms can be obtained for 
different quantities. Numerical simulations and analyt- 
ical estimations will also evaluate the effects of imper- 
fections in the quantum computer on the estimation of 
these quantities. 



II. HARPER AND KICKED HARPER MODELS 

The Harper model was introduced in 1955 to de- 
scribe the motion of electrons in a two-dimensional lattice 
in presence of a magnetic field. Its Hamiltonian reads 



HoiI,9) = cos(/) + cos(6l) 



(1) 
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FIG. 1: Phase space of the classical kicked Harper model: 
K — L ^ (Harper model) (upper left), K — L = 0.5 
(upper right), K — L — 1.5 (lower left), K = L = 2.5 (lower 
right) (10000 iterations of 256 classical orbits). One cell of 
size 27r x 27r is shown, the phase space being periodic. 

This Hamiltonian has been the subject of many studies 
(see for example 0, 0, 0, 0, 0), but its dynamics 
is somewhat restricted by the fact that it describes an 
integrable system. A generalization of this model has 
been introduced some time ago; it is called the kicked 
Harper model: 

H{I,e,t) = Lcos{I)+K cos{e)^ 5 {t-m), (2) 

m 

where m runs through all integers values and K, L are 
constants. This Hamiltonian reduces to ^ in the limit 
K = L 0, but has a more complex dynamics depend- 
ing on the parameters. Its dynamics between two kicks 
can be integrated to yield the map: 

r / = I + Ksm9 , . 

[0 = 6 - LsinJ ^ ' 

As in the case of the kicked rotator, there is a classical 
periodicity in both 6 and /. Thus the phase space is com- 
posed of cells of size 27r x 27r where the same structures 
repeat themselves. 

This map lj3Jl has been related to the motion of elec- 
trons in a perpendicular magnetic and electric fields, and 
also to the problem of stochastic heating of a plasma in 
a magnetic field. 

The quantization of yields a periodic Hamiltonian 
which after integration over one period yields a unitary 




FIG. 2: (Color online) Example of stochastic web in the 
kicked Harper model. Here K = L = 0.5, phase space is 
8x8 cells of size 27r x 2n, the figure shows positions after t = 
1000 iterations of 10® classical trajectories initially distributed 
according to a gaussian centered half a cell above the center 
with standard deviation ^/2tt/2^ ^ 0.0004. Color (grayness) 
shows density of points, from red (gray) (maximal value) to 
blue (black) (minimal value). 



evolution operator acting on the wave function ip 

{p^Uii^e-'^ cos(ftri)/r»g-«K cos(e)/h^^ ^^^^ 

where n = -iQd/dO and i){e + 2Qtt) = ^{9). 

This system has been the subject of many studies in the 
past few years, which focused on localization properties 
Jm ,2^ ,22, ^ ^ 21] , tunneling properties [H 
l29f . etc... 

In the limit K = L ^ Q the system is classically in- 
tegrable. For small iC, L, chaos begins to appear around 
separatrices, and spreads over larger and larger phase 
space areas as K,L increase (see Fig^. In the regime 
of small L, classical transport from cell to cell is pos- 
sible only in the very small chaotic zones around sepa- 
ratrices. For K = L, this network of thin chaotic zones 
surrounding large islands is called "stochastic web" (see 
Fig|2Jl . For intermediate values oi K,L, the phase space is 
mixed, with integrable islands separated by large chaotic 
zones. For larger K, L, classical chaos is present in most 
of the phase space (cf Fig^, and a typical trajectory will 
diffuse through the system. The quantum dynamics is re- 
lated to these classical properties, but shows some strik- 
ing differences. In the limit if = L — > 0, the system is 
integrable, wave functions are concentrated around clas- 
sical tori, but complexity manifests itself in the spectrum 
of the Hamiltonian, which is fractal ("Hofstadter butter- 
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FIG. 3: Map of delocalization in the {K, L) plane. Grayness 
represents the Inverse Participation Ratio ^ = 1/T,„\ijj{n)\'^ 
(IPR) , a measure of delocalization of states, from ^ = 1 (state 
localized on one momentum state) to ^ = Nh (totally delocal- 
ized state) (Ai'_ff is the dimension of the Hilbert space). Con- 
tour lines correspond to values of ^ ranging from 32 to 192 
by increments of 32, Nh = 2^, h/2n = (13 - VE)/S2 (actual 
value is the nearest fraction with denominator 2®). White 
corresponds to lowest values, black to maximal values of ^. 
Each ^ value is obtained by averaging over all eigenstates of 
the evolution operator f/ of Q|. 



fly" ) . For small K, L the motion of a quantum wave 
packet is dominated by the presence of classical invari- 
ant curves; the wave packet can spread in between these 
curves, or cross them by quantum tunneling. For larger 
K, L, in the regime of classical diffusion, as in the kicked 
rotator, a phenomenon similar to Anderson localization 
of electrons in disordered solids takes place. Through 
this phenomenon, called dynamical localization, a wave 
packet started at some value of momentum n will first 
spread, but contrary to classical trajectories this spread- 
ing will saturate. This corresponds to the fact that eigen- 
functions ijjain) of f7 in 10} in momentum space (they 
are called Floquet eigenfunctions since they correspond 
to the action of the evolution operator during one pe- 
riod) are exponentially localized. Their envelopes obey 
the law Tpain) ~ exp(— |n — m\/l)/Vl where m marks the 
center of the eigenstate and I is the localization length. 
This phenomenon is especially visible for moderate values 
of K, where all eigenfunctions are localized. For larger 
values of K, the system undergoes a transition: some 
eigenfunctions are still localized, but more and more are 
delocalized (ergodic) and spread over the whole system. 
This coexistence of localized and delocalized states gives 
rise to specific physical properties. Indeed, it is very dif- 



ferent from what happens in the kicked rotator model, 
where usually all states are localized once classical chaos 
is present (see for example or in the Anderson tran- 
sition (investigated in 30]) where the transition sepa- 
rates a regime where all eigenstates are localized from 
a regime where all are delocalized. In this regime of 
partial delocalization, an initial wave packet will spread, 
but a certain fraction of the total probability will re- 
main localized. In addition, the diffusion of probability 
in momentum space has been shown numerically to be 
anomalouSjWith an exponent depending on the parame- 
ter values [13, 122. 1231 . These properties are summarized 
by the phase diagram of Fig|3| Different quantities can 
be obtained in these different regimes with the help of a 
quantum computer. 

The phase space can be decomposed in cells of size 
27r X 27r. Its global topology depends on boundary con- 
ditions. For a system of size Nh, if the phase space is 
closed with periodic boundary conditions, with respec- 
tively Q and P cells in the 9 and n directions, then 
h = 2'kPQ /Nh- Therefore if one wants to keep h 
constant, the product PQ should be chosen such that 
PQ/Nh is the closest rational to fi./(27r). For most 
of the results of this paper, the phase space will be 
a cylinder closed in the 9 direction ((5 = 1) and ex- 
tended in the direction of momentum, and transport 
properties will be studied in the momentum direction, 
as in the kicked rotator. In this case h/ (27r) was set to 
l/(6 + l/((V5-l)/2)) = (13-\/5)/82 as in ^ to avoid 
unwanted arithmetical effects. The choice of a constant 
h implies that changing the number of qubits leads to 
increasing the size of phase space (number of cells) in 
the n direction. Only for the study of the stochastic web 
present at small K — L (subsection IV A) will the phase 
space be extended in both directions and its size (num- 
ber of cells) fixed. In this case increasing the number of 
qubits leads to smaller and smaller h. 



III. SIMULATING THE TIME EVOLUTION: 
THREE POSSIBLE ALGORITHMS 

The evolution operator @ is composed of two trans- 
formations which are diagonal in respectively the mo- 
mentum and position representations. This form is gen- 
eral for a family of kicked maps such as the kicked ro- 
tator, sawtooth map, and others. On a classical com- 
puter, the fastest way to implement such an evolution 
operator on a wave function of Nh components is to use 
the Fast Fourier Transform algorithm to shift back and 
forth between the n and 9 representations, and to imple- 
ment each operator by direct multiplication in the basis 
where it is diagonal. In this way, 0{NH^ogNH) clas- 
sical operations are needed to implement on a Nh- 
dimensional vector. On a quantum computer, it is pos- 
sible to use the Quantum Fourier Transform (QFT) to 
shift between momentum and position representations, 
using 0{{log Nh)^) quantum gates. In each representa- 
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tion, one has then to implement the multiphcation by a 
phase, e-^-^'=°''('^")/'i and e'"'^ . 

In the following we will envision three different strate- 
gies to implement these diagonal operators: 

• exact computation using additional registers to 
hold the values of the cosines 

• decomposition into a sequence of simpler operators 
which are good approximations during short time- 
slices 

• direct computation, the cosine function being ap- 
proximated by a (Chebyshev) polynomial 

The first one is in principle exact, but requires extra 
registers, and was already proposed in |0. The second 
one has some similarities with the one explained in [s^ 
for another system. The third one was not used in the 
context of quantum algorithms to the best of our knowl- 
edge, although the method is well known in computer sci- 
ence (see for example [sj for a recent use of this method 
to simulate many-spin systems on a classical computer) . 
We note that an approximate algorithm to simulate the 
kicked Harper for long time was used in ^32j : however, 
in that paper the aim of the authors was different, since 
they only wanted to construct efficiently a good approx- 
imation of the ground state wave function in order to 
use it for generating phase space distributions of other 
systems, and it is not clear that the method works for 
other purposes. We also note that the simulation of the 
Harper model on optical lattices was envisioned in [ssf . 
In the following discussions, we note by Uq the total num- 
ber of qubits including ancilla and workspace qubits, and 
N = 2"9 is the total dimension of the Hilbert space of the 
quantum computer. We note Ur with < riq the num- 
ber of qubits describing the Hilbert space of the kicked 
Harper model (i.e. the wave function evolved through Q 
is A^//-dimensional with Nh — 2"''), and Ug is the num- 
ber of elementary quantum gates used for one iteration 
of the quantum map I^J. 

A. Exact algorithm 

This approach is similar to the one taken in for the 
quantum simulation of the kicked rotator. In each repre- 
sentation, the value of the cosines is built on a separate 
register, and then transferred to the phase of the wave 
function by appropriate gates. 

If one starts with a A^jj-dimensional wave function 
— J2i2o~'^ CLi\(^i) in the 6 representation, with Nh = 
2"'', then the first step is to perform: 

Nh-1 Nh-1 

J2 a^\e,)\o) J2 a^\e,)\cos9,) 

1=0 i=0 

To this aim, the 2nr values cos(27r/2'') and sin(27r/2^), 
for j = l,..,nr are first precomputed classically with 



precision 2""'' with for example Up = 2nr using a re- 
cursive method based on Moivre's formula; then since 
9i = X]j=i (}ij2n/2^ with (3ij = or 1, one has: 

) = J|exp(zA,27r/2^) 

= Y[{cos{(3ij2n/23) + i sin(^y27r/2^)) 

This enables to compute | cosOi) for each 0i in n,. mul- 
tiplications by exp(i27r/2-') conditioned by the values of 
(3ij, needing in total O(n^) quantum gates. 

Then once the binary decomposition of cos 9i is present 
on the second register, conditional application of the 

one-qubit gates ^ ^ exp( iK2^^ /K) ) y^*^^*^^ ^^'^ state: 
aicyi]i{— iK cos{9 i)/h)\9i)\ cos 9 i) 

i=0 

Then the cosines in the last register are erased 
by running backward the sequence of gates that con- 
structed them, and one ends up with the state 
J2i2o~^ cos{9i)/h)\9i)\0), which is the result 

of the action of the unitary operator exp(—iK cos{9) / H) 
on 

Then the use of the QFT which needs 0{n'^) quantum 
gates shifts the wave function to the momentum represen- 
tation, and exactly the same technique as above enables 
to implement the operator exp{—iLcos{hh)/h) in 0{n^) 
quantum gates. A second QFT enables to go back to the 
9 representation. 

The whole process implements one iteration of the evo- 
lution operator U in 0{n'^) operations, with exponential 
precision. This algorithm is therefore efficient, and pre- 
cision can be increased exponentially at a cost of poly- 
nomial number of operations. On the other hand, the 
drawback of this approach is the need of several extra 
registers (one holding the values of the cosines, plus oth- 
ers for the workspace of the computation) and a relatively 
large number of gates. In the present status of experi- 
mental implementations of quantum computers, both the 
number of qubits and the number of gates that can be 
applied are very expensive resources. In the following, 
we will therefore expose two alternative strategies to im- 
plement U , which arc much more economical in the use 
of resources, but involve additional approximations. 



B. Slice method 

This technique enables to compute the operator U of 
Q without explicitly calculating the cosines. It approx- 
imates C/ by a sequence of many operators, each of them 
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being easier to compute. It can be viewed as "slicing" 
the operator into elementary operators. 

As above, we start with a iV//-dimensional wave func- 
tion \tp) = X^ilfo ^ ^A^i) ill the representation, with 
Nn = 2"''. In general, suppose we want to perform the 
operator 

Uk — e^'^^ (p ^) 

In the 9 representation, this operator is diagonal, so we 
just have to multiply each state by the phase e~''^ . 
First, we write 6 as 



— y 



(5) 



where the c?i's are the binary expansion of 6 and Nh 
2"''. If p = 2°m with m odd, then 



ptl = 



27rm 



/n,. — a— 1 \ 

da'^'' mod 2n 
1=0 / 



Thus Uk is equivalent to applying 

g — z/i: cos (m 6) 

on the rir ~ a first qubits. In the following, we will 
suppose that p is odd (a = 0) for the sake of simplicity. 

With the help of one ancilla qubit, let us perform the 
following sequence, where all gates are applied to the 
ancilla (initially set to |0)), except for Cfj which is the 
operator U applied on the principal register, controlled 
by the ancilla (the gate sequence is also displayed on 
FigEj): 

M{a,U) = HCuHe''^'^' HCu-^He''^"' HCuH 
This product is equal to 

M{a,U) = cos^ - - sm^ 



+1 sma 



2 2 2 



1 + ia 



2 



-Gy — I sm — - 



2 2i 
cr^ + 0{o?) for a < 1 



If we take U = e^^^ 

Af(a, [/) = \^'ia. cos {pd)a,_ +0(a^) 
since the ancilla qubit is in the |0) state, 

The kick operator can then be performed by ris 3> 1 
applications of Af(a, C/) 

Uk ~ Af(a, Uf" with a = — 



A more accurate expansion can be obtained by sym- 
metrizing Af (a, J7) 

= 1 + ^.-yl ^ ) +0(a3) 



Thus Uk ~ Af (a, U) up to order 2 in a. 

In this way, once a certain precision has been fixed, 
can be chosen such that the error is small enough. 

If we apply this strategy to the kicked Harper 
model, the method is therefore to first compute 
exp(— li^T cos(6')/?i) through the technique above (fc = 
K/h, p = \), then use a QFT to shift to the momen- 
tum representation. In the n representation, the opera- 
tor exTp{—iLcos{hn)/h) can be cast in the form above for 
h — 2T:'m/NH, with p ^ k ^ L/h and 9 — s- 2Trn/NH- 
The use of a QFT then shifts back the wave function to 
the 9 representation. 

The evolution of a A'jj-dimensional wave function with 
Njj = 2"'' through one time-slice is efficient, costing 
0{logNH) quantum operations. Indeed, for Us slices, one 
diagonal operator in is implemented in 4: + 2{nr^a) + 
{ns — l)(7+2{nr — a)) elementary gates, with a < rir. The 
number of slices fixes the precision. If one requires a fixed 
precision, independent of the number of qubits, then the 
whole method is efficient, iterating U in 0((log NhY) op- 
erations (the most costly operation asymptotically being 
the QFT). However, if one requires the precision to in- 
crease with Nh, then the method becomes less efficient. 
This algorithm is quite economical in qubits, since to 
simulate a wave function on a Hilbert space of dimension 
2""^, only Hq = rir + \ qubits are needed. One should 
note that for large number of slices, their computation 
dominate the computation time although asymptotically 
the QFT dominates. In all numerical simulations we per- 
formed, the slice contribution was indeed dominant. 

To precise the accuracy of the method, we show ex- 
amples of the localization length in the localized regime 
as a function of number of gates in FigEl The conver- 
gence with increasing number of slices (gates) is clearly 
seen, although for small number of gates oscillations are 
present. Data from = 7, 8 and rir = 9, 10 are close 
to each other due to the structure of the algorithm : in- 
deed, h/{2-K) is approximated by its closer approximants 
m/NHi and incrementing rir by one changes every other 
time the value of h. No major modification is seen in the 
numerical data for increasing rir, indicating that in this 
regime ris does not need to be drastically changed with 
rir- 

One may think that the spectrum is a much more sen- 
sitive quantity than the localization length. In FigEl we 
display the convergence for the spectrum of U for K and 
L small, in a parameter regime close to the fractal "but- 
terfly" visible for the unkicked Harper model (see Fig l24|l . 
The quantities displayed correspond to eigenphases Ea 
where ?7|V'a) = Q'^v{iEa)\^a) for some liAa)- The matrix 
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FIG. 4: Gate sequence for slices algorithm. Rz are Z rotations of angle —a. 
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FIG. 5: Localization length computed with the slice method 
over exact localization length Iq as a function of number of 
slices per iteration. Localization length is extracted after t — 
1000 iterations. Initial state is \tpo) = |0), with K = 1, L = 5, 
h/2-K = (13 — \/5)/82 (actual value is the nearest fraction with 
denominator 2""" with = — 1). 
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FIG. 6: Eigenphases of as a function of number of gates 
with the slice method: only 16 values are shown. Here = 6 
{uq ^nr + l),h = 27^/2^ K ^L = 10"^ 



C. Chebyshev polynomials 



of the operator L/ of Q is built by evolving through the 
slice method explained above the basis vectors, and then 
diagonalized. Convergence can be achieved with a few 
hundred time-slices. Due to numerical limitations, we 
cannot present data for different values of ri^, but we do 
not expect any drastic modification. 

In the subsequent sections, numerical simulations of 
this algorithm in presence of errors will be performed. 
To keep the computation time reasonable, we chose to 
use the slice method with 2 x 40 slices per iteration for 
transport properties (section IV) . Although the localiza- 
tion length is not exactly the correct one, the system is 
still localized and enables to study the variation of trans- 
port properties in presence of errors and imperfections. 
For computation of the spectrum (section V), we used 
2 X 100 slices per iteration. 



In this approach, one uses the QFT as in the preceding 
methods to shift back and forth between 9 and n repre- 
sentations. In each representation, the relevant operator 
is implemented by using a polynomial approximation of 
the cosines. Since polynomials can be implemented di- 
rectly through controlled operations, this avoids the use 
of additional registers. A commonly used polynomial ap- 
proximation rests on Chebyshev polynomials. 

Chebyshev polynomials (see for example [s^) are de- 
fined by the recurrence relation 



Ti{x) 

Tn{x) 



1 

X 



2xT„-i{x) — Tn-2{x) for n > 2 



They are bounded by —1 and 1 on [—1,1], with their 
extrema smoothly distributed over this interval. If f{x) 



is an arbitrary function on [—1.1], and we define for j 
0,..., M- 1 



A/-1 



cos 



fe=0 

Then, for large M 




cos 



is a very good approximation of f{x) on [—1, 1] 
If we truncate this formula to order m: 
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then the error is bounded by X]j=jri+i l*^*^! ^^"^ smoothly 
spread over [—1,1]. Practically, the c/c's are always 
rapidly decreasing, so the error term is dominated by 
|cm+i| and we can choose a small m while still keeping 
a good polynomial approximation of f{x). 

Let P{x) be a Chebyshev polynomial approximation 
of cos {tt{x + 1)). If one wants to perform the operator 

Uk = e^^^ '■^ on a A'^jj-dimensional vector with Nh = 
2"'' as in the preceding subsection, then for p = 1 



FIG. 7: Circuit for multiplying the quantum register 9 (sim- 
ple lines) by an odd classical number m (double lines) 



Uk can be decomposed as a product of operators of the 
form Ar{l3) = e^^^~^\ 
From lO, 



n 



f/fc only on the Ur ^ a first qubits. We then multiply 
the register by m before applying the cosine kick. 
Since m is co-prime with the dimension of the Hilbert 
space Nh = 2"'', this operation is unitary and can be 
performed without any additional qubit (for example 
with the circuit in Fig.|7l). 

If we choose a Chebyshev polynomial approximation of 
degree d, then the complexity of the algorithm is 0{nj.^^). 
This method is economical in qubits, and the precision 
of the approximation is easy to control. On the other 
hand, the complexity increases with the precision, and 
this can become prohibitive for very precise simulations. 
It is nevertheless quite efficient for fixed precision compu- 
tations, as can be inferred from the fact that it is actually 
the method used in classical computers to evaluate func- 
tions. 

In our numerical simulations, we found that a Cheby- 
shev polynomial of degree 6 was enough to get a very 
good approximation of the wave function. This demands 
a much larger number of gates than the slice method, 
^h-oM) the multi-controlled phase gate, which apply and scales badly with n,, in (here n, = rir since there 
the phase exp (z0) conditionally on the control qubits are no ancilla or workspace qubit). However, some of 
ji . . . jr (if an index is redundant, then it is counted only the control-phase gates have very small phases and are 
once), physically irrelevant. We can then choose a precision 

threshold and simply drop all the gates with phases be- 
low this threshold. The distribution of the phases of the 
gates computing the Chebyshev approximant of degree 6 
is displayed in the inset of Fig|H| 

This method of approximation is investigated in FiglSj- 
|31 The localization length as a function of number of 
gates is displayed in FigO for the same system parame- 
ters as in Fig 13 In Fig El we display the convergence for 
the spectrum, in the same regime as in FigEl 

In both cases, the convergence is good for maximal 
number of gates, showing that the polynomial of degree 
6 is indeed a good enough approximation in this regime 



Since the cZ^'s are binary digits, dj^ ---dj^ is equal to 
zero unless all terms are equal to one. If we denote by 



( ^\ Vl + -+3r 



Since all these gates commute, and since all the gates 
used for the construction of Ar are also present in the 
development of Ar' for r' > r, then all the terms of the 
polynomial P can be applied at the same time as the 
term of highest order by merging similar gates. 

If p 7^ 1, then p is split into p = 2°m with m odd, as 
in Ijlll B|) . The even part 2° is dealt with by applying 
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FIG. 8: Localization length computed with the Chebyshev 
method over exact localization length lo as a function of num- 
ber of gates. System parameters are the same as in Fig|3 with 
rir — 7 (dashed line), rir = 8 (dotted line), rir = 9 (full line) 
{uq = rir). Dashed horizontal line is / = Zq. Chebyshev poly- 
nomial of degree 6 is taken, keeping gates with the largest 
phases. Inset: number of gates as a function of their phase. 
Logarithms are decimal. 



of parameters. A good accuracy is achieved for a lower 
number of gates, implying that dropping the gates with 
the smallest phases can be an effective way to shorten 
the computation keeping a reasonable accuracy. Still, 
the data presented lead to the conclusion that even with 
the elimination of a large number of gates the method is 
clearly costlier in running time than the slice method to 
achieve similar precision. 



IV. TRANSPORT PROPERTIES: 
MEASUREMENT AND IMPERFECTION 
EFFECTS 

The three methods exposed above enable to simulate 
efficiently the effects of the evolution operator U of the 
kicked Harper model on a wave function. This produces 
the wave function at a given time. An important question 
concerns which quantities can be obtained through quan- 
tum measurement of the registers, and if the whole pro- 
cess including measurement is more efficient than classi- 
cal computation. A separate question but also related to 
practical efficiency of these algorithms is their stability 
with respect to errors and imperfections while running 
them on a realistic quantum computer. 

In this section, we will focus on the transport proper- 
ties of the wave function. We recall that for the kicked 
Harper model, for small K, L diffusion can only takes 
place on the small chaotic layer of the stochastic web. 
Then for larger K, L there is a regime of parameters 
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FIG. 9: Eigenphases of Q a function of number of gates 
with the Chebyshev method. System is the same as in Fig|S] 
Chebyshev polynomial of degree 6 is taken, keeping gates with 
the largest phases. An overall phase factor (global motion of 
eigenvalues) has been eliminated. 



where all eigenstates are localized, and another regime 
where localized and delocalized eigenstates coexist (see 
Fig|3J). In these different parameter regimes, we will show 
that quantities measuring localization properties and dif- 
fusion can be obtained on a quantum computer more ef- 
ficiently than on a classical device, although the gain is 
usually polynomial. We will then test the resilience to 
errors of these quantities obtained through the quantum 
algorithms, in particular through large-scale numerical 
computations. The error model we chose corresponds to 
static internal imperfections. Indeed, physical realization 
of a quantum computer will never be perfect, and some 
amount of disorder will always be present. In particular, 
interactions between qubits, which are needed to build 
the two-qubit gates, cannot in general be totally elimi- 
nated when they are not needed. These static imperfec- 
tions are not linked to interaction with the outside world; 
they have been shown to give important effects, which 
can be larger than the effects of noisy gates JJ, 35, SGj]. 
To model such errors, between each gate we require that 
the system evolves through the Hamiltonian 



(6) 



where the ct^ are the Pauli matrices for the qubit i and 
the second sum runs over nearest-neighbor qubit pairs 
on a circular chain. The energy spacing between the two 
states of a qubit is represented by its average value Aq 
plus a detuning 5i randomly and uniformly distributed in 
the interval [—5/2,5/2]. The detuning parameter 5 gives 
the width of the distribution near the average value Aq 
and may vary from to Aq. The couplings Ji repre- 
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sent the residual static interaction between qubits and is 
chosen randomly and uniformly distributed in the inter- 
val [—J/2, J/2]. We make the approximation that this 
Hamiltonian ^ acts during a time Tg between each gate 
which is taken as instantaneous. Throughout the pa- 
per, we take in general one single rescaled parameter 
e which describes the amplitude of these static errors, 
with e = STg = jTg. To probe the transport proper- 
ties of the kicked Harper model on a quantum computer, 
we chose to set H constant; in this way, changing the 
number of qubits is equivalent to changing the size of 
phase space (adding one qubit doubles the size of the 
phase space). The only exception is in the first follow- 
ing subsection (ncar-integrable regime), where the phase 
space volume is constant and h varies with the number of 
qubits. Throughout this section, effects of imperfections 
will be assessed using the slice method to implement Q . 
Therefore the presence of one ancilla qubit implies that 
n„ = n,. -|- 1 in all of this section. 




A. Near-integrable regime: stochastic web 

For K, L very small, the classical system is near- 
integrable : quantum transport is dominated by the pres- 
ence of invariant curves. Motion from cell to cell can take 
place only by tunneling effect, or by moving in the small 
chaotic zone around separatrices. In the case K = L, this 
small layer forms a "stochastic web" (see FigEJ which ex- 
tends in both 9 and n directions. A wave packet started 
in this region will slowly diffuse along this web. This 
process is best seen using quantum phase space distribu- 
tions, which allow direct comparisons between classical 
distributions such as the ones in Fig lll2l and quantum 
wave functions. 

The Wigner function [s^, HI] is an example of such 
quantum phase space distribution. However, it can take 
negative values, and only a smoothing over cells of area H 
gives non negative values. The use of a gaussian smooth- 
ing leads to the Husimi distribution (see e.g. js^) which 
in our case is defined by the formula: 



hie,n) 



2P 



n+NH/2 ^ 

1^(171) e 

ra=n-NH /2+1 



(7) 



where the gaussian for simplicity is truncated for val- 
ues larger than Nh/'2, ip{m) is the wave function in 
momentum representation, P (resp. Q) is the number 
of cells in the momentum (resp. position) direction, 
Nh = 2"'' is the dimension of the Hilbcrt space, and 
Q = NhO /{2TrQ). We note that methods to compute 
phase space distributions on a quantum computer were 
discussed in [11^2, 40]. 

In Fig^lwe show the spreading of a wave packet along 
the stochastic web for different numbers of qubits and 




FIG. 10: (Color online) Example of Husimi distribution of a 
wave packet spreading on the stochastic web; here K = L = 
0.5, h = 2iv X 64/2"'' (8x8 cells), initial state is a gaussian 
wave packet of area ft started half a cell above the center of 
the figure, after 100 iterations using 2 x 40 slices per iteration. 
Left: e = and from top to bottom Ur = 14, rir — 11, 
rir = 8 {nq = Ur + 1); right: rir — 14 and from top to bottom 
£ = 10~^, £ — 10~^, e = 10"*. Color/grayness is related 
to amplitude of the Husimi function, from zero (blue/black) 
to maximal value (red/white). Compare with the classical 
diffusion in FigEl 



different strengths of imperfections. In this picture, the 
size of the classical phase space is fixed, and the num- 
ber of qubits gives the value of h. A diffusion process 
is observed, which can be related both to the classical 
diffusion on the stochastic web (Fig|21) and to the effect 
of quantum tunneling from cells to cells. The diffusion 
constant is seen from Fig^|to depend on h; it also de- 
pends on K, L (data not shown) and is clearly different 
from the classical diffusion constant (compare the differ- 
ent times in Fig [21 and Fig llO|l . In this near-integrable 
regime, the tunneling process is quite complicated and 
was recently studied in p9j ]. In the same figure, one can 
see that with moderate levels of imperfections the exact 
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Husimi distribution is well reproduced by the algorithm. 

To probe transport properties in this regime, one can 
start a wave packet in the stochastic web and let it evolve. 
After a certain number of time steps, the diffusion con- 
stant can be obtained from measurement of the wave 
function. As the number of components of the wave 
function or of the Husimi distribution becomes expo- 
nentially large as rir increases, the best way is to use 
coarse-grained measurements : measuring only the first 
qubits adds up the amplitudes squared of many neigh- 
boring components and limits the number of measure- 
ments to a fixed value. This can be done to the wave 
function directly in the momentum or position represen- 
tation, or to the Husimi function provided all the val- 
ues are kept on a quantum register. For example, the 
Husimi-like function developed in |iT3 | can be obtained 
by modified Fourier transform from the wave-function, 
and allows the use of coarse-grained measurements. If 
one starts a wave packet on the stochastic web, it will 
diffuse according to the law (s(i)^) « Dgt, with s be- 
ing a distance in phase space and Ds the diffusion con- 
stant. Performing time evolution up to a time t* requires 
t* quantum operations multiplied by logarithmic factors. 
At this stage, a fixed number of coarse-grained measure- 
ments is enough to give an approximation of Ds- On a 
classical computer, one can truncate the Hilbert space up 
to the maximal dimension effectively used in the calcu- 
lation, which is of the order ^/t*. Propagating the wave 
packet will cost t* Vi* classical operations, after which Ds 
can be obtained. Therefore the quantum computation 
is polynomially faster than the classical one. Methods 
which use an ancilla qubit to measure the value of phase 
space distributions at a given point such as the ones in 
23 will necessitate extra measurements since they 
cannot be used to perform coarse-grained measurements 
efficiently. Still, by reducing if, L as rir is increased, one 
can keep the number of large components of the Husimi 
hmction of the wave packet of order Nh (instead of Nfj). 
In this case, the Husimi function measured on the ancilla 
qubit of (3^ is efficiently measurable. This is formally 
an exponential gain over direct classical simulation since 
measuring one component of the Husimi distribution at 
a fixed time t will be logarithmic in Nh- The same hap- 
pens for coarse-grained measurements at fixed t. Still, 
as h goes to exponentially small values the dynamics for 
fixed t will become very close to the classical one, so it is 
unclear which new information can be gained this way. 

To clarify the stability of these algorithms with respect 
to errors, in Fig^] we show quantitatively the effects 
of imperfections on the Husimi distributions for a wave 
packet spreading on the stochastic web for various num- 
bers of qubits and imperfection strengths. We computed 
the time scale th for various parameter values, th being 
the time (number of iterations) for which the error on the 
Husimi functions is half the mean value of that function 
on the stochastic web. 

The numerical data suggest the law 

th - C^/(e"<) (8) 
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FIG. 11: Effects of imperfections on the Husimi distribution 
of a wave packet spreading on the stochastic web ; here K — 
L — 0.5, h — 27rx 64/2"'" (8x8 cells), initial state is a gaussian 
wave packet of area h started half a cell above the center of 
the Fig|21 iterations are made by the slice method using 2 x 40 
slices per iteration. Straight line is the law (0 with a = 1 
and l3 = 1.23. Crosses corresponds to various values of e 
(10"'' < e < 10"") and Ur (5 < < 14, with Uq = rir + 1), 
averages were made over all Husimi components inside the 
stochastic web and up to 100 realizations of disorder for each 
e value. Inset : average relative error of the Husimi function 
Sh = — hd\) / {ho) on the stochastic web for e = 10"'* 
(dashed line), e = lO""'^ (dotted line), e = 10"^ (full line), 
rir = 10 {riq = rir -I- 1). Average is taken over all Husimi 
components inside the stochastic web and 10 realizations of 
disorder. Logarithms are decimal. 

with a = 1.02 ± 0.02 (compatible with a = 1) and 
13 = 1.23 ±0.09 with Ch ~ 0.007. This law is polynomial 
in both £ and n^, which indicates that even though indi- 
vidual values of the Husimi function can be exponentially 
small, the effect of imperfections remains small compared 
to these individual values for a polynomial time. This 
means that such quantities can be reliably obtained in 
presence of moderate levels of imperfections. More work 
is needed to understand the precise origin of the law ((SJ. 
We note that in fill where random noise in the quan- 
tum gates were used as main source of errors a similarly 
polynomial (but different) law was found for the relative 
error on the Husimi function. 



B. Localized regime 

When K is large enough for the chaotic zone to take 
most of the classical phase space, a classical particle will 
propagate diffusely in phase space. In contrast, for mod- 
erate values of the parameter K , all the eigenstates of 
the evolution operator t/ of l@J are localized (see Fig|21l. 
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FIG. 12: (Color online) Example of wave function in the 
localized regime, here K = 1, L ^ 5, h/{2-K) = (13 - ^)/82 
(actual value is the nearest fraction with denominator 2"'"), 
initial state is |V'o) = |0), after 1000 iterations using 2 x 40 
slices per iteration, = 8 (n, = + 1), from bottom to 
top e = (black, solid line), e = 10~^ (red, dashed line), 
e = 10"'' (green, solid line). In the center, the first two curves 
are superposed and indistinguishable. Logarithm is decimal. 



This localization is a purely quantum phenomenon due to 
interference effects and similar to the Anderson localiza- 
tion of electrons in solids. In this parameter regime, an 
initial wave packet will have projections on only a small 
number of exponentially localized eigenstates. Thus af- 
ter a few iterations of the map, the wave packet will stop 
spreading and stay in a region of momentum space of size 
given by the localization length. An example of such a 
wave function is shown in Fig ll2l 

In this regime, it is possible to measure the localization 
length I efficiently. Indeed, most of the probability is con- 
centrated in a domain of size If one performs a coarse 
grained measurement of the wave function, i.e. only the 
most significant qubits are measured, the number of mea- 
surements will set the precision in units of I. Thus once 
the desired relative precision is fixed, the number of mea- 
surements is independent of / or n^. Nevertheless, if one 
starts from an easily prepared initial wave packet, for ex- 
ample on a single momentum state, one has to evolve it 
long enough to reach a saturation regime where the wave 
hmction is spread on a domain of size « /. Classically, 
in the parameter regime where the system is chaotic, the 
dynamics is diffusive w Dt with a diffusion con- 

stant D which depends on parameters. One can expect 
the wave packet to follow for short time this diffusive be- 
havior which will stop when a spreading comparable to 
the localization length is reached. In this case, the wave 
packet needs to be evolved until a time t* sa P /D. Clas- 
sically, one needs to evolve a vector of dimension ~ I until 
the time t*\ this needs ^ classical operations. On a 
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FIG. 13: Example of IPR with imperfections as a function 
of time, in the localized regime. Parameters values are the 
same as in Fig ll2l = 8 {uq = Ur + 1), from bottom to top 
£ = , £ = 10"'^, e = 10"'', £ = 10"^. Data from £ = and 
£ = 10"^ are indistinguishable. 

quantum computer, once the precision is set the three al- 
gorithms above need only a logarithmic number of gates 
to perform one iteration, so the total number of gates is 
^ . This gives a polynomial improvement for the quan- 
tum algorithm. It is known that in the delocalized phase, 
the wave packet can spread ballistically for some regimes 
of parameter. If this extends to short times and to the 
localized regime, then the gain becomes quadratic. 

In Fig ll2l an example of a localized wave function is 
shown for different imperfection strengths. At e = 0, the 
exponential localization is clearly visible, the exponential 
decay being leveled off at very small values (~ 10"'^°) 
only by numerical roundoff. For larger values of e, the 
localized peak is still visible with the correct amplitude, 
but a larger and larger background is visible, until the 
peak disappears. 

To analyze in a more precise way the effects of imper- 
fections, we have to specify the observable that is used 
to get the localization length. On a classical computer, 
different data analysis can be used to calculate the local- 
ization length from knowledge of the wave function. A 
first way consists in extracting the second moment of the 
wave function ((An)^), which gives an estimate of I once 
the saturation regime is reached. One can also compute 
the inverse participation ratio (IPR) C = 1/S„|V'(«)|^. 
For a wave function uniformly spread over M states this 
quantity is equal to M, and therefore it also gives an es- 
timate of the localization length. At last, I can be mea- 
sured directly by fitting an exponential function around 
maximal values of ip. 

For an exact wave function, all three quantities give 
similar results. On a quantum computer, they may 
have very different behavior with respect to imperfection 
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FIG. 14: IPR as a function of imperfection strength in the 
localized regime. Parameters values are the same as in Fig |12l 
with rir = 7 (full line), rir = 8 (dotted line), n,. — 9 (dashed 
line), Tir = 10 (long dashed line) (n, — n^ + l). Averages were 
made over up to 10 realizations of disorder. Inset: Fidelity 
as a function of imperfection strength in the localized regime, 
with same parameter values and line codes as in the main 
figure. Logarithms are decimal. 



strength. Indeed, it was shown in general [lll|41| that the 
second moment is exponentially sensitive to the number 
of qubits in presence of imperfections, making it a poor 
way to get information about transport properties. The 
IPR was shown to be polynomially sensitive to both 
number of qubits and imperfection strength. Still, the 
IPR may be difficult to measure directly on a quantum 
computer. On the other hand, the direct measurement 
of I by fitting an exponential curve on a coar se-g rained 
measure of the wave function was shown in '4l' to be 
an effective way to extract / from a quantum computa- 
tion of the wave function. It is therefore interesting to 
study the behavior of both latter quantities with respect 
to imperfections. 

In Fig^l the time evolution of the IPR is shown for 
different values of the imperfection strength. For e = 0, 
the wave packet first spreads for t < t* then the IPR 
becomes approximately constant and close to the local- 
ization length. For larger values of e, the wave packet 
spreads to much larger parts of phase space, but the IPR 
still saturates after some time to a value which depends 
on £ and Uq. 

The average value of this saturation value is shown in 
Fig^las a function of e for different values of Uq. FigllSI 
shows the localization length obtained from curve-fitting 
for the same wave functions. For large enough values of 
e, the IPR grows very quickly, in a manner which seems 
exponentially dependent on Uq. The result of the curve- 
fitting strategy is roughly similar, but shows an interme- 
diate regime {e ~ 10^'' for our data) where it is still rea- 



FIG. 15: Localization length as a function of imperfection 
strength in the localized regime. Parameters values are the 
same as m Figini with averages made over up to 10 realiza- 
tions of disorder. Logarithm is decimal. 



sonably close to the exact value while the IPR is already 
quite far off. This can be understood qualitatively from 
the data shown on Fig^| Indeed, the effect of moderate 
static imperfections is to create a larger and larger back- 
ground over which the localization peak is superimposed. 
The IPR is sensitive to the presence of this background, 
while by its very definition the curve-fitting strategy iso- 
lates the localization peak from the background and is 
therefore more robust. The data presented in Figll5lshow 
that this peak keeps its shape with relatively good accu- 
racy until its final disappearance, even though a large 
chunk of its amplitude has been transfered elsewhere by 
imperfections. The inset of Fig^] shows the fidelity of 
the same wave functions, (/(t) = |(-0(i)|V'e(i))P where 
IV'(i)) is the exact wave function and |V'e(i)) the one in 
presence of imperfections). It is interesting to note that 
the localization length and IPR can be quite well repro- 
duced even for values of e where the fidelity is already 
quite low. 

A more precise analysis can be developed from the 
effect of imperfections on the eigenstates of the unper- 
turbed evolution operator [/ in Q. These eigenstates 
\ipa) can be written as a sum over momentum states |m), 
which coincide with quantum register states of the quan- 
tum computer when the system is in momentum repre- 
sentation : 



N„ 

E 



(9) 



In the localized regime, the eigenstates l-^a) are lo- 
calized with localization length I, therefore the are 
significant only for ~ I values of to, with average value 



13 



Using perturbation theory, one can estimate the 
typical matrix element of the imperfection Hamiltonian 
© between eigenstates. For the first term of ©, this 
gives: 



Vt 



typ 



{'^pb\'^S^cr^Tgng\^pa) 

Nh 



i=l 



(10) 



where Nh - 
which U acts, 



: 2""^ is the dimension of Hilbert space on 
Tg is the time for one gate, and the term 
due to Ao in © is not taken into account since it can 
be eliminated easily. This estimate (jlOf) is an approxima- 
tion, since the action of JB)) is separated from the action 
of U and in reality they are intertwined and do not com- 
mute. In (|10|1 . only ~ / neighboring quantum register 
states are coupled through Uq terms of different detun- 
ing Si (with random sign). This term therefore gives on 
average eng^Jn^j \fl. The second term of © in the same 
approximation will be the sum of Uq terms, each coupling 
one state |to) with another state differing by two neigh- 
boring qubits |n) = |r7i + r). So a state j-^a) is coupled 
significantly only to states localized at a distance r 
in momentum from |'0a)- Therefore the same estimate 
applies, and overall one can estimate Vty-p ^ sn^ 



-q/Vl. 

One can suppose that the IPR will become large when 
perturbation theory breaks down. This happens when 
Vtyp is comparable to the distance between directly cou- 
pled states Ac. From the arguments above, one expects 
that one state is coupled to ~ / states so that this dis- 
tance is Ac ~ l/l. The threshold when IPR or localiza- 
tion length become large is therefore Vtyp ~ Ac which 
corresponds to: 



Ci/(ng^/n^\fl) 



(11) 



where Ci is a numerical constant and Ug is number of 
gates per iteration, Uq number of qubits, / the localiza- 
tion length. Fig ll6l is compatible with this scaling, with 
C\l\fl ~ 0.3. We note that this threshold is similar to 
the threshold for the transition to quantum chaos pre- 
sented for a quantum computer not running an algorithm 
in 

When perturbation theory breaks down, it is usu- 
ally expected from earlier works on quantum many-body 
physics [35I 0j that the system enters a Breit-Wigner 
regime where the local density of states is a Lorentzian 
of half- width F w 27r|Vfypp/Ac according to the Fermi 
golden rule. This implies that the IPR grows like F/A„ ~ 
e'^n^nqN, where A„ ~ ~ l/N is the mean level 

spacing {N — 2Nh since there is an ancilla qubit). This 
is not confirmed by the data shown on Fig ll7l which sug- 
gest that the IPR scales like e. This indicates that in our 
system we are in a regime different from the golden-rule 
(Breit-Wigner) regime. 
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FIG. 16: Critical value of e (error strength) as a function of 
parameters for K = 2, L = 27, with other parameter values 
the same as in Fig |12l Sc is defined by a saturation value of 
IPR twice the unperturbed value. Averages were made with 
up to 10 realizations of disorder. Solid line is the formula 
lllll . Logarithms are decimal. 
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FIG. 17: IPR as a function of e. Parameters values are the 
same as in Fig |16l Solid lines correspond to the dependence 
^ oc e . Logarithms are decimal. 



Such a regime is present for large perturbation strength 
in many-body systems. It is indeed known that for large 
enough values of the couplings, the system leaves the 
golden rule regime and enters a new regime where the 
local density of states is a gaussian of width given by 
the variance a. The variance can be approximated by 
cr^ ^ y^(,j„ V,?,,, e^riina. In this regime the IPR is 



given by cr/A,, 



-g'-q- 



n-qN, 



which is consistent with 
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data from Fig^] and FigEI This regime is known to 
supersede the golden rule regime for T > a, which should 
therefore be the case for our system. This implies that 
the relevant time scale for the system to remain close 
to the exact one is 1/a. For the largest values of Ug, 
the data in Fig 1 171 show some departure from this law, 
which may be due simply to statistical fluctuations (the 
averaging is made over more instances for smaller Uq), or 
a shift toward the golden rule regime for large Uq. 

The scaling laws obtained in this regime show that for 
e < Ec, with Ec given by (|ll|l . the system is still localized 
in presence of imperfections, and the localization length 
is close to the exact one. In this case, the localization 
length is correct for very long times, much longer than for 
example the fidelity decay time. For larger e, the system 
with imperfections is delocalized. We still expect it to be 
close to the exact one up to a time ^ 1/a ^ 1/ {eng,Jfi^). 
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C. Partially delocalized regime 

For larger values of if at L fixed, the system enters 
a partially delocalized region. In this regime, there is a 
coexistence of localized and delocalized eigenstates. An 
initial wave packet will have significant projections on 
all delocalized eigenstates but only on neighboring local- 
ized eigenstates. After a certain number of time steps 
(kicks) the part corresponding to delocalized states will 
spread in all the system, while the localized part will 
remain close to the initial position. Fig 1 181 shows an ex- 
ample of a wave packet initially at n = after 100 it- 
erations in this regime, displaying an exponential peak 
corresponding to localization superimposed on a plateau 
which spreads with time to larger and larger momentum. 
It is known that the spreading of the wave packet in this 
regime (for large enough time) is ballistic away from the 
line K = L and diffusive on this line. 

In this regime, as above a coarse grained measurement 
can give the localized part with moderate accuracy, thus 
enabling to compute the localization length. As in the 
preceding part, the gain over classical computation will 
be polynomial. As concerns the delocalized part of the 
wave function, it seems at first sight that getting infor- 
mation on it is difficult, since it takes very long time 
to reach its saturation distribution (it has to spread dif- 
fusively or ballistically through the whole system), and 
this distribution itself is spread over the exponentially 
large system. Still, after a time large enough for the 
wave packet to spread beyond the localization length, the 
structure of the wave function can be seen very clearly 
from coarse-grained measurements whose number is on 
the order of the localization length. Once such coarse 
grained measurement has been performed, and the lo- 
calization length found by fitting an exponential func- 
tion around the maximum, the relative importance of 
the plateau can be found by subtracting the localized 
part. Even though the plateau has not yet reached its 
final distribution, its integrated probability is related to 



FIG. 18: (Color online) Example of wave function in the 
partially delocalized regime. Here K = 2, L = ^, h/{2Ti) = 
(13 — v'5)/82 (actual value is the nearest fraction with de- 
nominator 2"""), initial state is I'i/'o) = |0), after 100 iterations 
using 2 X 40 slices per iteration, = 8 (uq = -I- 1), from 
bottom to top e = (black, solid line), e = 10"*^ (red, dashed 
line), £ — 10~^ (green, solid line). In the center, the first 
two curves are superposed and indistinguishable. Logarithm 
is decimal. 



the number of eigenstates which are delocalized at these 
parameter values. This information enables to monitor 
the transition precisely for different values of K and L, 
a non trivial information as seen from Fig|31 The num- 
ber of operations for classical and quantum algorithms 
are the same as for the localization length, and there- 
fore the same polynomial gain can be expected. Another 
quantity which can be readily obtained is the quantum 
diffusion constant. Indeed, away from the line K = L, 
it is known that a quantum wave packet initially local- 
ized in momentum will spread anomalously (ballistically) 
with the law (n^(i)) « Dat^- Classically, estimating the 
diffusion constant requires to simulate the system un- 
til some time t* . This requires to evolve ~ quantum 
states until the time i*, making the total number of op- 
eration ~ (t*)^- On a quantum computer, one time step 
requires a logarithmic number of operations, so the to- 
tal number of operations is ~ t* [t* iterations followed 
by a constant number of coarse-grained measurements), a 
quadratic gain compared to the classical algorithm. Close 
to the line K = L, the quantum diffusion becomes nor- 
mal with the law {■n?{t)) « Dnt. In this regime, the 
same computation gives a number of operation ~ (i*)'^/^ 
classically compared to ~ t* for the quantum algorithm, 
with still a polynomial gain. Such computations can give 
quite interesting results, in particular to specify precisely 
which kind of diffusion is present in the vicinity of the 
line K — L, a question which is not definitively settled. 

Effects of different strengths of imperfections can be 
seen in Fig^l For moderate values of e, a flat back- 
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FIG. 19: Example of IPR in presence of imperfections as a 
function of time in the transition regime. Here K — 4, L — 5, 
h/{2Tv) = (13 — ^/5)/82 (actual value is the nearest fraction 
with denominator 2"""), initial state is \^o) = |0) with 2 x 40 
slices per iteration, rir — 8 {uq — rir + 1), from bottom to top 
e = 0, £ = 10"'', e = 10"**, e = 10"^. Data from e = and 
£ = 10~^ are indistinguishable. 



FIG. 20: Critical value of e (error strength) as a function of 
parameters for K — 10, L = 27 with other parameter values 
the same as in Fig llSI Sc is defined by a saturation value of 
IPR twice the unperturbed value. Averages were made with 
up to 10 realizations of disorder. Solid line is the formula 
11211 . Logarithms are decimal. 



ground of larger and larger amplitude is created by the 
imperfections. When this background reaches the val- 
ues of the plateau due to dclocalized states, information 
on these delocalized states is lost, but the localized peak 
remains until e is large enough to destroy it. This is 
visible also on Fig 1 191 which displays the time evolution 
of a wave function in this transition region. The data 
for e = show the spreading of the wave packet due 
to delocalized eigenstates; the IPR does not reach the 
dimension of Hilbert space since part of the amplitude 
does not spread due to localized states. For intermediate 
values of e, the spreading concerns more and more of the 
total amplitude, increasing the IPR, until a large enough 
value of e is reached and the wave function is completely 
delocalized. 

In this regime, the analysis of the preceding section 
should be modified. Indeed, a certain fraction /3 of the 
Floquet eigenstates lipa) of U (unperturbed) in Q are 
not localized. For these delocalized states, the c™ of ^ 
have small nonzero values ~ for all m. The esti- 

mation Vtyp ~ eng^Jn^l Vl for the typical matrix element 
of the imperfection Hamiltonian © between eigenstates 
I ■(/'a) and I ■(/'ft) remains correct only if |?/;a) a-nd It/;^) are 
both localized. 

If It/iq) and I'f/'fc) are both delocalized, one has c™ ^ 
Cfe ~ l/V^-ff in H10|l for most m, n. This implies that 
the quantities X]m=i '^o'c™: previously of order 1/%/! be- 
comes ~ l/V-ZV// (sum of terms of order ~ l/A^// 
with random signs). This modifies the estimate for Vtyp'- 
with the same reasoning as in the localized case, one has 



If one of the states \il}a) and is localized and the 
other one delocalized, then X]m=i '-o'c™ i^ the sum of I 
terms of order ^ ^/{\/1\/Nh) with random signs, which 
is of order ^ 1 / Nh ■ This gives the same estimate 
Vtyp ^ ^ng^Jn^l \JN}i for the matrix element as if both 
states were delocalized. 

Therefore if a proportion /3 of the unperturbed Flo- 
quet eigenstates are delocalized, both localized and de- 
localized eigenstates will have matrix elements of order 
Vtyp ~ eng^n^l \fNu with (iNu other eigenstates. This 
will be the dominant effect, since these couplings lead 
perturbation theory to break down much earlier than for 
the purely localized system. Indeed, Vtyp is compara- 
ble to the distance between directly coupled states ~ 
XjNn ~ 1/iV (since N = 2Nh) for eug^/^/N - 1/iV, 
which corresponds to: 

Sc ~ C2/{ng^VN) (12) 

where C2 is a numerical constant, rig is the number of 
gates per iteration, Ug the number of qubits and N — 
2"' the dimension of the Hilbert space of the quantum 
computer. Fig l20l is compatible with this scaling, with 
C2 ~ 7.4. 

In this regime, the critical interaction strength drops 
therefore exponentially with the number of qubits, in 
sharp contrast with the localized regime. This effect has 
been noted for a different system in |43| . and is similar 
to the enhancement of weak interaction in heavy nuclei 
|44j . The physical mechanism is that the much smaller 
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FIG. 21: IPR in presence of imperfections as a function of K 
in the transition regime. Here L = 27, fi/ {2tv) = (13 — \/5)/82 
(actual value is the nearest fraction with denominator 2"""), 
initial state is \tpo} = |0), IPR is shown after 100 iterations 
using 2 X 40 slices per iteration, n,- — 8 (uq — rir + 1), with 
averages made over 10 realizations of disorder. 
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FIG. 22: Shift of the transition point due to imperfections 



as a function of imperfection strength and 



Parameters 



values are the same as in Fie l21l with averages made over up 
to 10 realizations of disorder. Solid line corresponds to the 
dependence AK^ oc eng^/n^N . Logarithms are decimal. 



coupling term between states is compensated by the even 
smaller distance in energy between coupled states. This 
result implies that even for moderate number of qubits, 
a small interaction strength is enough to modify enor- 
mously the long-time behavior of the system: saturation 
values of the IPR are very much affected by the pertur- 
bation, much more so than in the localized regime. How- 
ever, for short time the behavior of the system should be 
close to the unperturbed one, implying that the measures 
suggested to get interesting information, such as relative 
size of the plateau and diffusion constants can still be 
accessible. 

FigE] shows examples of the growth of the IPR as 
a function of K and imperfection strength. In the par- 
tially delocalized zone, the figure shows a growth of the 
IPR with K , which is strongly affected by imperfections. 
An interesting quantity is the value of the transition 
point between localized and delocalized states. In sys - 
tems such as the Anderson model investigated in |30l| . 
the transition point is well-defined, since all states are 
localized or delocalized on one side of the transition. 
In the case of the kicked Harper model there is some 
arbitrariness in the definition. We chose as transition 
point the value Kc (at L fixed) for which the IPR is 
Nh/4: — N/8 (even for a totally delocalized state, the 
IPR is actually often Nh /2 = iV/4 instead of Nh due 
to fiuctuations). In the partially delocalized regime, the 
IPR at fixed K should grow with e. If the system is 
in the Breit-Wigner (golden-rule) regime, IPR should 
grow as r/A„ where A„ ~ 1/A^ is the mean level spac- 
ing and r « 27r|Vtj,pp/Ac ~ e^n^n^. We therefore ex- 



pect the transition point to move with imperfections as 
r/A„ ~ e'^n'^nqN. On the contrary, in the gaussian 
regime, the IPR grows like (t/A„, where a eng^Jfi^. In 
this case the transition point should move as eng^Jn^N . 

Figim shows the data numerically obtained for the 
shift of the transition point due to imperfections. It 
indicates that Aifc ~ sng^JrT^N agrees with the data, 
whereas e'^n'^nqN is a much less reasonable scaling vari- 
able (data not shown). The data therefore seem to in- 
dicate that in the partially delocalized regime as in the 
localized regime, the IPR grows as eugy/n^N, as does the 
shift of the transition point. This result is in sharp con- 
trast with the findings of "30] for the Anderson transition, 
which was shown to scale polynomially with the number 
of qubits. In our case, the presence of delocalized state 
coexisting with localized states makes the delocalization 
much easier in presence of imperfections. 

FigESI shows the scaling of the IPR as a function of 
e. For small values of Uq, the IPR without imperfections 
is already a large fraction of Hilbert space dimension, so 
data are meaningful only for > 9. Still, data shown 
on Fig|531seem to indicate that the regime where ^ cx e 
is present, confirming that the system is in a gaussian 
regime rather than in the golden rule regime. 

The scaling laws obtained in this regime show that 
there is an exponentially small value Sc given by l|12|) 
above which imperfections destroy the localization prop- 
erties of the system. In particular, the transition point 
is exponentially sensitive to the number of qubits. This 
sharp difference between localized and delocalized regime 
can be easily seen on experiments: the long time behav- 
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FIG. 23: IPR as a function of e for A' = 10, L = 27 with other 
parameter values the same as in Fie llSI Solid lines correspond 
to the dependence ^ oc e . Logarithms are decimal. 



lor of the system will be very different in both cases. 
Still, the algorithms presented can be useful in delocal- 
ized regime in presence of imperfections, even for e > Ec- 
Indeed, the system should remain close to the exact one 
up to a time ~ 1/cr ^ 1/ {eng^Jfi^) as in the localized 
regime, so measurability of physical quantities will even- 
tually rest on the comparison of this time scale with the 
time for the system to show the delocalization plateau. 
On the contrary, in the localized regime for moderate lev- 
els of imperfections the localization length can be mea- 
sured for very long times. 



V. SPECTRUM: MEASUREMENT AND 
IMPERFECTION EFFECTS 

Another type of physical properties which can be ob- 
tained through quantum simulation of the kicked Harper 
model concerns the spectrum of the evolution operator 
U . This spectrum has been the focus of many studies (see 
s-g-H^I)- it shows multifractal properties, and transport 
properties (localized or delocalized states) are reflected in 
the eigenvalues, as well as dynamical properties (chaotic 
or integrable states). Additionally, for small K — L, 
this spectrum will be close to the famous spectrum of 
the Harper model ( "Hofstadter butterfly" ) , which shows 
fractal properties [ij], as can be seen in Figl^ 

To get information about eigenvalues, we can use the 
phase estimation algorithm. This algorithm, at the heart 
of the Shor algorithm, proceeds by transforming the state 
Y,t l^)IV'o) into J2t l^il^^Vo)- Then a QFT of the first 
register will give peaks at the values of the eigenphases 
of U. To be efficient, this process should be applied to 
operators U for which exponentially large iterates can be 
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FIG. 24: Eigenphases of the Harper operator @ as a function 
of ft for = L = 10~^ Ur = 8. 



obtained in polynomial number of operations. In |45| it 
was suggested that one can obtain approximate eigenval- 
ues exponentially fast provided one starts from an initial 
state iV'o) already close to an eigenvector. In the case 
at hand, we do not know how to get exponentially large 
iterates in polynomial time, nor how to build a good ap- 
proximation of the eigenvectors without knowing them. 
We therefore suggest a third strategy, which is more gen- 
erally applicable than the ones above, but does not yield 
an exponential gain. 

We first build the state J^t^o^^ l*)IV'o)j where \ipo) is 
an arbitrary quantum state on a Hilbert space of dimen- 
sion Nh = 2"'' which can be efficiently built, for example 
it can be the state 2""''/^ \n), which can be obtained 
from |0) with iir Hadamard gates. Once the state |0)|?/'o) 
is realized, it can be transformed with Hadamard gates 
on the first register into 2^"''/^ ^^'J,""^ l*)IV'o)- We have 
seen that the evolution operator U can be implemented in 
poly(logA'^ff ) operations by the three strategies exposed 
in section III. Therefore we can apply powers of U on 
the second register controlled by the value of the first 
register. This yields 2-''-/^J2t l^)|C^'V'o) in 0{Nh) op- 
erations, up to logarithmic factors. A QFT of the first 
register will yield peaks centered at eigenvalues of the 
operator U. Thus measurement of the first register will 
give an eigenvalue of U with good probability in 0{Nh) 
operations including measurement. A drawback of this 
approach is that peaks have additional probabilities on 
nearby locations, and since the number of eigenvalues is 
Nl{^ measuring the precise shape of all peaks will be in- 
efficient {0{Nfj) operations). A more precise, although 
slower method is to use amplitude amplification (a 
method derived from the Grover algorithm Q) to zoom 
on a small part of the spectrum. This enables to get the 
precise values of all eigenvalues in a given interval. The 
total cost will be 0{Nh\/Nh) operations. This methods 
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FIG. 25: Eigenphases of the evolution operator U of ^ as 
a function of imperfection strength. The slice method is used 
with 2 X 100 slices to compute the operator. The 16 eigen- 
phases closest to are shown. Here Ur = 6 (uq — + 1), 
h — 2tv/2^ (actual value is the nearest fraction with denom- 
inator 2^), K — L = 10"''. An overall phase factor (global 
motion of eigenvalues) has been eliminated. Logarithm is dec- 
imal. 



which uses Grover's search on phase estimation can be 
seen as a process reverse to quantum counting (where 
phase estimation is used on the Grover operator). 

Calculating the spectrum by direct diagonalization of 
a Nh X Nh matrix such as the one of the operator U of 
Q requires in general of the order of Nfj classical op- 
erations. However, in the case of the operator [/ of Q 
there is a faster classical method similar to the quantum 
phase estimation algorithm: one iteration of U can be 
computed classically in 0{Nh) operations (up to loga- 
rithmic factors) by using the classical FFT algorithm to 
shift between n and 9 representations, and multiplying by 
the relevant phase in each representation. Iterating this 
process Nh times and keeping each intermediate wave 
function costs 0{N^) operations. Then a FFT enables 
to get the spectrum of U with 0{Nh log A^h) operations. 
This method was advocated in By] for getting the spec- 
trum of the kicked Harper model. Therefore it is possible 
to get the spectrum classically in 0{Nfj) operations up 
to logarithmic factors. Thus the quantum algorithms ex- 
plained above {0{Nh) operations for one eigenvalue with 
unknown precision, 0{Nh\/Nh) for all eigenvalues in a 
given small interval) realize a polynomial gain compared 
to the classical ones. It is important to note that al- 
though the number of operations needed is only polyno- 
mially better in the quantum case, the spatial resources 
are exponentially smaller (logarithmic number of qubits 
compared to the number of classical bits). 

The Fig l25l27l show the spectrum of the kicked Harper 







-1 



LU 



O 



-3 



-4 



- 














- 











-10 -9 -8 -7 

log(e) 

FIG. 26: Average error (in units of mean level spacing) of 
computed eigenphases through the slice method as a func- 
tion of imperfection strength; parameters are the same as 
in Fia' 1251 and averages were made over all eigenvalues and 
over 10 realizations of disorder. Dashed line corresponds to 
oc e. Logarithms are decimal. 



model in presence of errors for both slice and Chebyshev 
methods. The error model chosen is the static imper- 
fection Hamiltonian © as in the preceding section. The 
evolution operator was computed by evolving basis states 
in presence of errors and then diagonalizing the result- 
ing operator. The spectrum shown corresponds to small 
K = where the spectrum is close to the "Hofstadter 
butterfly" , as can be seen in Fig[21 Only 16 eigenvalues 
are shown. Overall phase shifts due to errors were elim- 
inated since it seems reasonable they can be estimated 
and compensated. It is clear from the data presented 
that eigenvalues are much more sensitive to strength of 
errors than transport properties. Numerical limitations 
prevented us to find the scaling in Uq of error effects, but 
Fig l26l281 show the scaling with respect to e at constant 

In the case of the slice method, the average error on 
the eigenvalue is clearly linear in e. We think this cor- 
responds probably to a perturbative regime, since small 
values of e are involved. For the Chebyshev method, our 
data indicate that a lower level of errors is needed than 
in the slice method to get good accuracy. This could 
have been expected, since we established in Section HI 
than this method necessitates more gates for a similar 
accuracy, and each gate introduces errors. The scaling of 
errors with respect to e indicates the law AE ^ with 
a « 1.3. 
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algorithm can be meaningful, but a careful choice of the 
measured quantities should be done. For the different 
quantities computed, the slice method was shown to be 
more efficient and resilient to errors than the Chebyshev 
method, although the latter is similar to the method used 
in classical computers to evaluate functions. 

Our results show that interesting quantum effects 
such that fractal- like spectrum, localization properties, 
anomalous diffusion are already visible with 7-8 qubits. 
We therefore believe that such algorithms could be used 
in experimental implementations in the near future. 
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FIG. 27: Eigenphases of the evolution operator U of ^ as 
a function of imperfection strength. The Chebyshev method 
is used; a Chebyshev polynomial of degree 6 is taken, keeping 
all gates. The 16 eigenphases closest to are shown. Here 
rir = 6 {uq — Tir), ft — 27r/2® (actual value is the nearest 
fraction with denominator 2®), K = L — 10~^. An overall 
phase factor (global motion of eigenvalues) has been elimi- 
nated. Logarithm is decimal. 



VI. CONCLUSION 
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In this paper, several quantum algorithms were pre- 
sented which enable to simulate the quantum kicked 
Harper model, a complex system with relevance to cer- 
tain physical problems. The comparison showed that 
while the slice method and the Chebyshev method are ap- 
proximate, they are much more economical in resources 
than the exact simulation. It was also shown that dif- 
ferent transport and spectral properties can be obtained 
more efficiently on a quantum computer than classically, 
although the gain is only polynomial. Numerical simula- 
tions enabled us to precise the effect of numerical errors 
on these algorithms, and also to evaluate the effects of 
imperfections. The results show that depending on the 
regime of parameters, the same quantity can be stable or 
exponentially sensitive to imperfections. In general, in 
presence of moderate amount of errors the results of the 



FIG. 28: Average error (in units of mean level spacing) of 
computed eigenphases through the Chebyshev method as a 
function of imperfection strength; parameters are the same 
as in Fig |27l and averages were made over all eigenvalues. 
Dashed line corresponds to A_E oc e^ '^. Logarithms are deci- 
mal. 
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